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Oscillations of rotating magnetised neutron stars with 
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ABSTRACT 

We investigate the oscillation spectrum of rotating Newtonian neutron stars endowed 
with purely toroidal magnetic fields, using a time evolution code to evolve linear 
perturbations in the Cowling approximation. The background star is generated by 
numerically solving the MHD equilibrium equations and may be nonspherical by virtue 
of both rotation and magnetic effects; hence our perturbations and background are 
fully consistent. Whilst the background field is purely toroidal, the perturbed field 
is mixed poloidal-toroidal. From Fourier analysis of the perturbations we are able 
to identify a number of magnetically-restored Alfven (or a-) modes. We show that 
in a rotating star pure inertial and a-modes are replaced by hybrid magneto-inertial 
modes, which reduce to a-modes in the nonrotating limit and inertial modes in the 
nonmagnetic limit. We show that the r-mode instability is suppressed by magnetic 
fields in sufficiently slowly rotating stars. In addition, we determine magnetic frequency 
shifts in the /-mode. We discuss the astrophysical relevance of our results, in particular 
for magnetar oscillations. 
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1 INTRODUCTION 

The identification of Soft Gamma Repeaters (SGRs) and Anomalous X-ray Pulsars (AXPs) as magnetars has reignited interest 
in the role played by magnetic fields in stellar dynamics. A class of neutron star (NS) with surface fields of ~ 10 15 gauss, 
magnetars are the most highly magnetised stars yet observed. The SGRs are of special interest since in addition to their regular 
gamma-ray flares, they are also susceptible to occasional giant flares over three orders of magnitude more energe tic than the 



regular flares. In the tails of the se flares, high-frequency quasi-periodic oscillations (QPOs) have been observed (jlsrael et al 
20051 : IWatts fc Strohmaverll2007l — these may be the first direct detections of neutron star oscillations. These QPOs should 



provide a valuable probe of the physics of the NS interior, but this requires an understanding of stellar oscillations in the 
presence of a strong magnetic field. 

Neutron stars are not the only stars where magnetic effects may be important. The influence of a star's magnetic field on 
its oscillation spectrum can be gauged from the ratio of its magnetic energy to the gravitational binding energy, M/|W|; this 
suggests three classes of star where one should take account of the star's magnetic field: in addition to NSs, there are also the 
rapidly-oscillating type- A peculiar (roAp) stars and magnetic white dwarfs (MWDs). 

The earliest studies of magnetic star oscil lations were driven by the discovery of ~ 10 4 gauss fiel ds — relatively strong 
for a main-sequence star — in some Ap stars ( Chandrasekhar fc Limber 19541 : Ledoux fc Simor] 1957 ). Later, some of these 



stars were found to be oscill ating at high frequency — the roAp stars — m o tivating a number of stud ies of magnetic effects 
on high-frequency p-modes (jUnno et al. I llQ8Ql : iDziembowski fc Goodelll996l : iRincon fc Rieutordl Eooi h In addition to roAp 
stars, some white dwarfs have strong (~ 10 9 gauss) magnetic fields; these ha ve shown no evidence of pulsation , perhaps due 
to magnetic suppression of the (/-modes observed in weaker-field white dwarfs) Wickramasinghe fc Ferrario 2000h . Finally, the 
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internal dynamics of neutron stars will be considerably affected by rotation as well as their strong magnetic fields, leading to 
interest in magnetic r-modes (jMorsink fc Rezanial 120021 ). 

Many publications to date have reported on analytic studies of magnetic stellar oscillations, necessitating considerable 
simplications to the problem: typically the model used is an incompressible star with a force-free background magnetic field. 
Some modern work on the problem has been inspired by the observation of magnetar QPOs, and this has tended to be numerical 
( Sotani et aHl2008l : Isotani fc Kokkotadl2009l : Icerda-Duran et al.ll2009l ). with the advantages that more sophisticated physics 
can be modelled (for example, compressible and relativistic stars). 

This paper is an investigation of the mode spectrum of magnetic stars through time evolution of the perturbed MHD 
equations. Here we consider only purely toroidal background fields, but in future we intend to investigate purely poloidal and 
mixed field stars too. We begin by perturbing the full system of equations to yield a set of equations for the background and 
another set for the perturbations. We improve on previous work on magnetic oscillations by solving the background equations 
in a self-consistent manner, so that the star is in exact equilibrium (up to some small numerical error), rather than simply 
using spherically symmetric density profiles and analytic results for the magnetic field configurations. We next discuss details 
of our time-evolution code and test its accuracy. After this, we present results for oscillations of stars with magnetic fields 
and rotation. Finally, we interpret these in the light of observed magnetic star oscillations and compare them with previous 
work on this subject. 



2 BACKGROUND AND PERTURBATION EQUATIONS 

We model a neutron star as a self-gravitating, rotating, magnetised polytropic fluid with perfect conductivity. The system is 
then governed by the equations of ideal magnetohydrodynamics (MHD): 

P I ~T" + ( v ' v ) v + 20 x v ] = -VP - pV$ - ptl x (O x r) + — (V x B) x B, (1) 

\Ot J 47T 

V 2 $ = 4vrGp, (2) 
t = -V-(pv), (3) 

^=Vx(vxB), (4) 

P = kp\ (5) 

together with the solenoidal constraint V-B = on the magnetic field. Here v denotes the part of the fluid's velocity field which 
is not rigid rotation CI; all other symbols have their usual meanings. Throughout this paper we work with 7 = 2 polytropes 
exclusively, as a simple approximation to a neutron star equation of state. We consider linear Eulerian perturbations of this 
system by making the standard ansatz that each physical quantity has a zeroth-order background piece and a first-order 
perturbed piece; e.g. the density is written as p — po + 5 p. 

We assume that our background star is stationary and rigidly rotating, so that f2 is zeroth-order and v first-order. 
Equations <[3j and Q become trivial and we are left with 

= -VP - p V$<, - p Q n x (O x r) + — (V x B„) x B , (6) 

47T 

V 2 $ = 4irGp , (7) 

Po = kpl (8) 

Making the additional assumption of axisymmetry one may show that this system of equations splits into two cases: one 
where the magnetic field is purely toroidal and a second mixed-fi e ld cas e (with pure-poloidal fields as a limiting case). Details 
of the solution of these equations are given in lLander fc Jones! (|2009h : we use the code described therein to generate the 



background configurations used here. Here we merely note that our background configurations are fully self-consistent, with 
rotation, magnetic fields and fluid effects in equilibrium. In contrast to other work on magnetic oscillations, our background 
star need not be spherical, but may be distorted by rotational or magnetic effects, or a combination thereof. 

For the perturbation equations, we work in the rotating frame of the background and write our equations in terms of 
the perturbed density Sp, the mass flux f = pov an d a magnetic variable (3 = po^B. We additionally make the Cowling 
approximation — neglecting the perturbed gravitational force — to avoid the computational expense of solving the perturbed 
Poisson equation. Our perturbations are then governed by seven equations: 

Po^r = -7JW5p-20xf + ( (2 ~ 7)7P ° Vp - -1(V x B ) x B ) 5 p +-L(VxB )x/3+-^(Vx/3)xB () -— L(Vp ( ,x/3)xB 

at V Po 47T / 47T 47T 47rpo 

(9) 
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dS P 
dt 



= -Vf, 



9/3 -Vx(f xB„)-^x(f xB ). 



(10) 
(11) 



dt p 

If we rewrit e these equa tions in terms of SP = -^-Sp and set the magnetic field to zero they reduce to the perturbation 
equations of Jones et al. 1 20021 ). 

We next decompose each perturbed quantity in the azimuthal angle <j>; for example Sp becomes: 



8p(t, r, 9, <f>) = 5p^(t, r, 6) cos m<j) + Sp m (t, r, 8) sin m<j). 



(12) 



This reduces our problem from a 3D to a 2D one, at the expense of doubling the number of equations: we now have evolution 



equations in fourteen variables: f^, f^ 1 ,Sp ± ,(3^ z ,(3^ 1 , (3^ for some value of m to be specified at the start of the evolution. 



2.1 Boundary conditions 

Rotational and magnetic forces will serve to distort the star's density distribution away from spherical symmetry and hence 
complicate the treatment of perturbations at the stellar surface. To avoid these complications we replace the radial coordinate 
r with one fitted to isopycnic surfaces, x = x(r,8); even a nonspherical surface will be defined by one value x = R. With the 
background density being a function of x alone, we have po(x = R) = and hence 

f (x = R) = (5{x = R) = 0. (13) 

Finally, the Lagrangian pressure perturbation AP is zero at the surface by definition. Relating this to the Eulerian perturbation 
we have 

6P + £-VP = at the surface. (14) 

Using JBJ, we see that VPo may be written as two terms proportional to po and a term involving the magnetic current 
V x B/4-7T. Both density and current are zero at the stellar surface and so VPo must also vanish there. This yields our last 
surface boundary condition: 

SP(x = R) = 0. (15) 

Our boundary conditions allow us to evolve the interior magnetic field perturbations of our star, but not oscillations 
of the exterior. By contrast, one would expect magnetic perturbations in a physical neutron star to reach the surface and 
produce electromagnetic radiation extending through the exterior. Whilst our treatment of the surface does not account 
for this, we believe that it is the most that can be done using the equations of perfect MHD: in an infinitely-conducting 
polytropic star, a magnetic field that extends to the surface has a corresponding Alfven speed ca = y B 2 /4irp which becomes 
superluminal as p — > 0, and infinite when p = (i.e. the stellar surface and exterior). Dealing with the surface and exterior 
thus requires extra physics: a stellar model more sophisticated than a polytropic fluid with perfectly electrical conductivity. 
One could employ a low-density numerical atmosphere for the exterior, or assume that the field is confined or matches to 
some simplified crust — but these are merely numerical conveniences rather than good models of actual NS physics. In reality, 
perfect MHD ceases to be a good approximation close to the surface of a NS, where resistive effects become important and 
the full equations of electromagnetism should be used. The stellar surface is not fluid but an elastic crust; and the exterior 
will have a magnetosphere region rather than a dilute, uniform 'atmosphere'. 

Needless to say, a credible model star which included all these effects would give an oscillation spectrum closer to that 
of a real neutron star than the one we study here. In lieu of such a model, however, we treat oscillations over the fluid, 
highly-conductive interior of the star only. With magnetic fields being strongest here and ~ 99% of the NS's mass consisting 
of a fluid interior, we anticipate that dynamics in this region would dominate the star's oscillation spectrum; and hence that 
our treatment is a reasonable first attempt to understand oscillations in real NSs. 

Next we look at the conditions at the centre of the star. Since we deal with m > perturbations in this study, we should 
enforce a zero-displacement condition: 

5P(x=0) = , f(x- = 0) = /3(x=0) = 0. (16) 

In general, a parity analysis of the perturbation equations of a fluid star shows that the variables may be classed 
according to their equatorial symmetry — either odd (the perturbation is zero along the equator) or even (its theta-derivative 
is zero). This division of the perturbations allows us to reduce our numerical domain to just one 2D quadrant and enforce the 
perturbation symmetry at the equator as another set of boundary conditions. 

Analysing the perturbation equations for the (unmagnetised) rotating fluid problem, one finds that the perturbation 
variables may be divided into the two symmetry classes {ff 1 , ft, $P ± } an d {fg 1 }- In the case of a background star with 
a pure poloidal field these classes are augmented by magnetic variables, viz. {/jr, ft , Sp ± , /3^}, {fg 1 , Pr, Pt}- Note that 
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although the background field is pure-poloidal, the perturbed field will still be mixed poloidal-toroidal. For a pure-toroidal 
background the magnetic perturbations are again mixed, but they fall into different symmetry classes from perturbations of a 
pure-poloidal star: , Sp ± , f3^r , j3^} and {f^tP^}- It follows that whilst we may separately treat perturbations on either 
a pure-poloidal or pure-toroidal background, the perturbations of a mixed-field background will have no definite equatorial 
symmetry. Investigating this latter group of perturbations requires an extended numerical domain consisting of an upper and 
lower quadrant. For this paper we concentrate only on oscillations of stars with purely toroidal background fields, postponing 
the pure-poloidal and mixed-field cases to future work. 



3 NUMERICS 

As described above, our numerical domain is one 2D quadrant of a circle, with x £ [0, 1] and 9 £ [0,n/2]; by symmetry and 
through a (^-decomposition this domain is sufficient to investigate behaviour over the whole 3D, potentially nonspherical, star. 
Upon decomposing in <f), equations Q, (|10p and become a system of fourteen perturbation equations, which we evolve 
using a MacCormack predictor-corrector algorithm. 

Our code is an extended version of the one described in Passamonti et al. 1 20091 ). where the authors demonstrated its 



accuracy and long-term stability for barotropic and stratified stars, in the nonmagnetic case. As in their work, we employ a 
fourth-order Kreiss-Oliger dissipation; this is an extra term added to the equations to damp spurious higher-order oscillations, 
i.e. those generated by the code's finite differencing. The magnitude of this term is resolution-dependent, so that it vanishes 
in the infinite-resolution continuum limit. 

In addition to this dissipation, two further tricks are required to ensure stability and accuracy of magnetic evolutions. To 
stabilise the numerical evolution of the magnetic field, we first note that if the electrical resistivity r] is non-zero, the induction 
equation gains an extra term: 

^=Vx(vxB)-t|Vx(VxB). (17) 

By including this second term (at a small magnitude) we are able to suppress instabilities which arise from evolving the 
magnetic field. As for the Kreiss-Oliger dissipation, this artificial resistivity is added in a resolution-dependent manner, so it 
reduces to the correct (i.e. physical) continuum limit. We find that a very small value of r] is sufficient to improve long-term 
stability, but has negligible physical effect on our evolutions, since it acts over a far longer timescale than any others in our 
problem. 



3.1 Divergence cleaning 

Finally, for the long-term accuracy of the code we need to ensure that the perturbed magnetic field remains solenoidal. This 
is guaranteed in the continuum limit, since the divergence of the induction equation is 

9( -^ B) = V ■ V x (v x B) = 0, (18) 

but in practice numerical error will be introduced from the finite grid resolution. It is important to 'clean' the field of this 
class of numerical error, since i t has been shown that a numerically-generated monopolar field gives rise to a spurious extra 
force |Brackbill fc Barnes|[l98r} ). 



There are various approaches to divergence cleaning for numerical schemes. A review of these may be found in lDedner et al 



( 2002h . where in addition a new constrained formulation of MHD is proposed, where the condition V ■ B = is coupled to 



the induction equation through an auxiliary function; we repeat their argument below. 

In the continuum limit the induction equation states that the vector <9 t B has a divergence- free part only, whereas a 
general vector can be decomposed into curl-free and divergence-free parts. Our discretised induction equation will no longer 
preserve this divergence-free aspect exactly and accordingly we add a curl-free term —Vip to the RHS, with tp being some 
unknown function. We then couple our augmented induction equation to a relation for ip: 

d t B = V x (v x B) - V?/> (19) 
V(i>) = -V-B (20) 

where T> is some linear differential operator. The Euler equation and the equation of mass conservation are unaffected. We 
now take the divergence of the first relation and the time derivative of the second: 

9t(V-B) = -VV (21) 
d t V(i>) = -9t(V-B) (22) 

which we combine to see that 

dtD(ip) = V 2 tp. (23) 
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The choice of D determines the way in which divergence errors are removed. The three basic types of cleaning are 
elliptic, para bolic and hype r bolic — so named because they entail solving a heat equation, wave equation or Poisson equation, 
respectively. iDedner et al. pioneer a mixed hyperbolic-parabolic approach, which they find to be superior to the 

simpler divergence-cleaning methods since it allows for errors to be propagated out of the star (hyperbolic cleaning) whilst 
simultaneously being damped (parabolic cleaning). The third method, elliptic cleaning, has the serious disadvantage that it 
requires the repeated solution of the (computationally expensive) Poisson equation; the mixed-hyperbolic scheme only adds 
the modest expense of having to evolve one more quantity — the function %p. 

Hyperbolic-parabolic divergence cleaning involves defining O by 



C h C P 



which leads to a telegraph (damped-wave) equation for ip: 



d tt iP = -^d t iP + c 2 h V 2 ip. 
Within the code, we implement this divergence-cleaning method through the evolution equation 



(24) 



(25) 



(26) 



together with our modified induction equation (|19|l . Following I Price fc Monaghanl <|2005l l we take ch, the divergence- wave 
propagation speed, to be related to the sound c s and Alfven ca speeds through the relation: 



Ch = yci + c 2 A . 

The other coefficient is physically the inverse of the decay timescale r: 

cl 1 



(27) 



(28) 



which iPrice fc Monaghanl ([20051 1 argue is not universal, but rather should be modified to suit some lengthscale A specific to 
the problem, i.e. 

where a is a dimensionless parameter. Using this result, we take A to be the radial grid spacing Ar in our code. Finally then, 
our evolution equation for the function ip is 



d t ip 



iP - (cj + ci)V • B. 



(30) 



To close the system we need to give appropriate boundary conditions and initial data. For the latter we simply set 
ip(t = 0) = — this is reasonable because the initial data is divergence-free and so the variable ip, associated with the 
monopole part of the magnetic field, should be zero initially. 

For the boundary condition at the surface, we choose the Sommerfeld outgoing wave condition on ip: 

d t Tp = -y/c i B +(? A (ft-V + VO- (31) 
This result is for a spherical surface, but it still gives satisfactory cleaning in the case where the background star is spheroidal. 



3.2 Testing the code 

Since we already have confidence in the performance of the code in the nonmagnetic limit (see Passamonti et al. ( 20091 ) 



for details), we now test its accuracy and convergence properties with the inclusion of magnetic effects. To this end, we 
wish to monitor the divergence of the magnetic field and the total energy of the system (which should be conserved in the 
continuum limit). Since the background is stationary i ts total energy is autom atically conserved; in addition the background 
field was constructed in a divergence- free manner (see lLander fc Jonesl (|200Sh .1. Therefore it suffices to check conservation of 



the perturbed energy and the value of V • <5B. 

To get a measure of the divergence of <5B over the whole star, rather than at individual points, we would like to evaluate 
V • <5B volume-integrated over the star. However, this quantity is first-order and so integrates to zero, so instead we define a 
'monopole energy' 

© = ^ j {V-SB) 2 dV (32) 

where R is the stellar radius; this energy should be insignificant in comparison with the perturbed magnetic energy SM. For 
each evolution used to generate results for this paper, the divergence of <5B was monitored through the dimensionless quantity 
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Figure 1. We determine the order of convergence of our code by evaluating the total perturbed energy 5£ over time; in the exact, 
continuum limit this quantity will not deviate from its initial value. The left-hand plot shows the deviation of 8<£ from its initial value 
for (r, 8) grids of 32 X 30 and 64 X 60 points. From these we confirm that the order of convergence Oconv of the code is equal to 2, as 
intended (see right-hand plot). O CO nv is only plotted for t 10, since at early times the numerical values of <5(E cross the continuum 
value, causing O CO nv to oscillate rapidly. The background configuration for these tests was a star with rotation rate Q/y/Gp = 0.238 and 
with an average magnetic field strength <B>= 2.87 X 10 16 G, evolved for 30 /-mode oscillations. 



D /SM. Typically this quantity was found to be of the order ~ 0.01, rising to ~ 0.1 in the case of very strong fields and rapid 
rotation rates. 

Next we use conservation of (perturbed) energy to test the order of convergence of our code, using the fact that in the 
limit of infinite resolution energy should be exactly conserved. The total energy <5£ of the perturbations is given by 

6<E = 5T + 6W + 5U + SM, (33) 

i.e. the sum of the perturbed kinetic T, gravitational W , internal U and magnetic M energies. Since we are making the 
Cowling approximation, SW = and we are left with 

S<E = ST + SU + SM = f (IpoW\ 2 + ^S p 2 + ^-\SB\ 2 ) dV. (34) 
J \2 j£o 8yr / 

This is in agreement with equation (C5) of Friedman fc Schutz ( 19781 ) in the case of adiabatic perturbations within the 
Cowling approximation, but with an additional magnetic energy term. Note that all these terms are second-order quantities; 
the first-order energy vanishes for a background in equilibrium, and in our case each piece of it (for example the magnetic 
energy term proportional to Bo • SB) is automatically zero by virtue of our (^-decomposition. 

To evaluate the convergence ratio, we monitor the evolution of the high-resolution energy <5fl;64x6o(£) and the medium- 
resolution energy <5fi:32x3o(i), comparing these with the initial value of the energy <5£(0). In the continuum limit 5<£ will have 
no time-dependence and will be equal to its initial value for all time. Hence we are able to use this exact result to define a 
convergence ratio 

Oconv = : log — }f — . (35) 

log 2 \6<£.64x6o(t) - 5<£(0)J 
In figure [T] we evaluate O CO nv over time, confirming that the code is second-order convergent. 



3.3 Nondimensionalising 

Throughout the code we employ variables which have been made dimensionless through division by a suitable combination of 
powers of gravitational constant G, central density p c and equatorial radius r eq . For example, a dimensionless mode frequency 
a is related to the physical one a (with units of rad s -1 ) through the relation a = a/^/Gp c ; the conversion is the same for 
rotational frequency £7. Since dimensionless frequencies of this form are common in oscillation mode literature we use these 
throughout this paper. Dimensionless magnetic field strengths, however, are less likely to be familiar and so we quote these 
in terms of gauss. 

When we use dimensional quantities they are for a neutron star with canonical parameters: an equatorial radius of 
10km (in the non-rotating, unmagnetised case) and a mass of 1AM§ (where M© is solar mass). The relationship between 
dimensionless frequencies a (equivalently Cl) and their physical counterparts is only weakly dependent on f2 and B — and 
hence is roughly linear, with 

£7 [Hz] w 1890a. (36) 
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Finally, we note that in our dimensionless units, the Keplerian (break-up) velocity Qk ~ 0.72. When we plot sequences of 
modes in rotating stars, we typically track the modes up to Q/Qk ~ 0.95; that is, rotation rates 95% of the break-up velocity. 



4 RESULTS 

4.1 Classes of oscillation mode 

Using spherical polar coordinates, a general perturbation may be decomposed with respect to the basis (Y; m e r , rVYj m , e r x 
VY; m ), where Yj m = Yi m (6, <fi) are the usual spherical harmonics. The first two of these terms transform by multiplication by 
(— 1)' under parity inversion r i— > — r, with the latter one transforming as (— . This enables us to classify modes based on 
their parity: those whose sign is given by (—1)' under parity interchange as termed polar modes, whilst those transforming as 
(— are called axial modes. Hybrid modes, consisting of a sum of axial and polar pieces, are termed axial-led or polar-led 
based on whether the lowest-Z (i.e. l = m ) term of the mode is axial or polar, respectively. 

For fixed m, lLockitch fc Friedman! (|l999h found that inertial modes are not characterised by a single I, but have an 



angular dependence consisting of a sum of spherical harmonics Yi m (9,<f>). However, in all cases they found there was some 
threshold value Iq, such that the amplitude of Yi m contributions for I > l was found to drop off rapidly. Following their work, 
we label modes using the notation where k distinguishes between different modes with the same lo- 

A non-rotating, unmagnetised fluid star has only one class of oscillation modes if the perturbations are assumed to have 
the same equation of state as the background star. These are the p-modes, whose restoring force is pressure fluctuations; the 
lowest-order p-mode (i.e. the one with a nodeless eigenfunction) in each series is termed the fundamental mode, or /-mode. 
The non-axisymmetric p-modes are degenerate in the absence of rotation and magnetic fields; each p-mode has the same 
frequency for fixed m. These modes are polar in nature. 

With a rotating background star, a Coriolis force term enters the equations governing the perturbations, which removes 
the m-degeneracy in the p-modes. The Coriolis term is the restoring force for a new class of modes: the inertial modes, which 
we term i-modes. In general i-modes are mixed axial and polar even in the slow-rotation limit, but one class of them are 
purely axial in this limit: the r-modes. With the barotropic equation of state we employ here, the only r-modes which exist 
are those with l = m. 

Finally, magnetic fields also induce a class of oscillation mode, restored by the Lorentz force. We term them the Alfven 
modes, or a-modes. In addition to generating a n ew class o f modes, the Lorentz force can lift degeneracies of nonradial 
oscillations, causing a splitting in mode frequencies 1 Cox 1980). The addition of the Lorentz force term in the Euler equation 



for the perturbations should produce shifts in the frequencies of the p, r and i modes from their unmagnetised values. 



4.2 Initial data 

Having performed a (^-decomposition of our perturbation variables, we fix a value for the azimuthal index m for each evolution. 
However, since we have no restrictions on ^-dependence, oscillations for a variety of I m will typically be excited for arbitrary 
initial data — hence we are able to study several modes at once. The nature of these modes is then determined from analysis 
of their eigenfunctions, as well as by tracking them back to a regime where we already know their properties. 

Whilst any initial data will excite a variety of modes, we choose different starting perturbations depending on whether we 
wish to investigate axial/axial- led or polar/polar-led modes. For polar modes we provide a density perturbation whose angular 
dependence is given by an ordinary spherical harmonic; for axial modes we use a 'magnetic' spherical harmonic-dependent 
velocity perturbation . In b oth cases the radial dependence is a Gaussian profile. More details about these choices may be 
found in I Jones et al. |2002[ ). 



4.3 Mode spectrum of a non-rotating magnetised star 

In this section we present results for nonrotating stars, since the mode spectrum is simpler, leaving rotating stars to the next 
section. We begin by investigating the new class of modes present with the addition of a magnetic field, the a-modes. Results 
are presented for both polar and axial a-modes. 

Let us begin by considering where in the frequency spectrum these modes could be expected. Now, any mode frequency 
will be proportional to some characteristic wave speed. For fluid modes like the /-mode, the frequency should be proportional 
to the sound speed c s ; similarly the a- mode frequencies should be proportional to the Alfven speed ca ■ Accordingly the ratio 
of frequencies should scale as 
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Figure 2. Typical FFT results for a pair of nonrotating stars, one magnetised and the other unmagnetised. We plot mode frequency 
a (in a dimensionless form) against PSD, the power spectral density. We see that the /-mode frequencies are very close in both cases. 
With no Coriolis force there are no inertial modes, therefore any peaks at lower frequency than the /-mode must be either noise or 
Alfvcn modes. We identify the lowest-frequency spike in the magnetic FFT as noise, since there is a corresponding unphysical peak in 
the nonmagnetic FFT. The following peaks in the magnetised-star FFT, however, have no analogue in the nonmagnetic FFT and so we 
identify these as Alfven modes. The duration of the evolution was sufficient to resolve around 100 Alfvcn oscillations. 



where the angle brackets represent a volume average. Now 



CA 




and so 



2Vtt7 <P> 



(38) 



(39) 



t c A I <B> 

We find from our background code that a nonrotating unmagnetised 7 = 2 polytrope with a mass of 1.4M© and radius R = 10 
km has a volume- averaged pressure <P> of 3.10 x 10 34 dyn cm" 2 . Using this value and <B>= 10 16 G to nondimensionalise, 
we find that 



^~90x 



1/2 



<B> 



,3.10 x 10 34 dyn cm" 2 ,/ V 1Q16 G ,' 
With the value of < P> varying little with magnetic field strength, we assume that it is a constant and that oj/<J a scales only 
with <B>. It then follows that we should expect a a to be roughly 1/90 of Of for a 10 16 G field, but 1/9 of 07 for a 10 17 G 
field. This part of the spectrum is dominated by inertial modes in the case of unmagnetised rotating stars, but in the absence 
of rotation we may be confident that any oscillations at lower frequency than the /-mode are associated with the magnetic 
field — see figure [2] 

Now, with a oc< ca > and ca = B/^/inp, it follows that a a oc B, provide d that magnetic chang es to the density 



Lander fc Jonesl ( 20091 )). To summarise 



distribution are higher order (which should be true for all but very high fields — see 
a-modes should scale linearly with field strength and appear as oscillations with lower frequency than the /-mode. With these 
expectations, we now turn to numerical results from our time-evolution code. 

In figure [3] we track a number of Alfven mode frequencies up to averaged-field strengths of order 10 17 gauss. For axial 
initial data and fixed m we find a single I = m mode, whilst polar initial data excites two lo = m + 1 modes for a given m. In 
all cases, we see that as expected there is a near-linear relationship betwee n a a and <B>. The iden tification of the a-modes 
is based on analysis of their eigenfunctions, using the numerical method of IStergioulas et al.l (|2004 ). The labelling used here 
anticipates the results of the next section, where we track these modes for increasing rotation rate. 

At the start of this section we showed that the a- mode frequency should vary linearly with <B>, and this appears to be 
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4.5 



shift in u f (%) 




<B> 2 [gauss 2 ] 

Figure 4. The shift in /-mode frequency due to magnetic effects (for nonrotating stars), for m = 2, 4, 6. On the y-axis we plot percentage 
increase in from its unmagnetised value; we see that this shift appears to depend quadratically on <B>. The apparent deviation from 
this dependence, visible in the weakest-field results, is attributable to numerical errors in these very small frequency shifts. 



borne out by our results. We now quantify this dependence and the deviation from it. By looking at the weak-field results from 
our code (where the relationship should be closest to linear), we determine the constants of proportionality in the relationship 



10 16 G 

finding that |c = 0.033, |ci = 0.030, |c 2 = 0.045, \c = 0.086, |ci = 0.069, \c* = 0.090, %c = 0.146, £ci = 0.127, \c 2 = 0.150. For 
all our results the linear relationship (|41[) . with the numerically-established constants Jjct, agrees to within 8% of the mode 
frequencies we have extracted from our evolutions — and in most cases the difference is less than 5%. 

Finally in this section, we look at the shift in the frequency of the fundamental mode upon the addition of a magnetic 
field to the star. This mode is restored by perturbations in the fluid pressure P in the unmagnetised case, so we anticipate that 
in the magnetic problem the restoring force is perturbations of total (fluid plus magnetic) pressure, P + B 2 /8ir. The magnetic 
shift in erf, then, should be proportional to B 2 — but since magnetic pressure is very modest in magnitude compared with 
fluid pressure, we expect the frequency shift to be small. For example, using our canonical model star, the magnetic pressure 
is ~ 1% of the fluid pressure at 10 17 G. We confirm these expectations in figure|4] In all cases 07 is increased by the inclusion 
of magnetic effects, but the shifts are only around a couple of percent even for <B>~ 10 17 . The relative shift appears to be 
more pronounced for higher-m oscillations. 
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Figure 5. Magnetic shift of the m = 2 /-mode frequency for rotating stars. Since the shift is very small we take a very highly magnetised 
background star, with <B>= 1.17 X 10 17 G for comparison with the nonmagnetic sequence of results. We find that as the rotation rate 
Q increases, magnetic effects become less significant. 




<B>= 1.44 x 10 16 G 
<B>= 2.87 x 10 16 G 
<B>= 5.71 x 10 16 G 




Figure 6. Illustrating the hybrid incrtial-magnetic nature of modes in a rotating magnetised star. When Q = there is a pure I = m = 2 
a-mode, which is split into co- and counter-rotating modes by the effect of rotation. The counter-rotating mode frequency approaches 
the nonmagnetic r-mode frequency as Q increases, while the corotating branch tends to zero frequency. The left-hand plot compares 
the 2 a mode with the r-mode, whilst the right-hand plot shows that the nature of the hybrid mode depends on the ratio M/T; when 
<B> is larger, the jd-mode frequency approaches the r-mode frequency more slowly. Modes are tracked up to Q ~ 0.7 in dimcnsionlcss 
units, which is over 95% of the break-up velocity. The irregular parts of the curves may correspond to avoided crossings with other 
magneto-inertial modes. 



4.4 Mode spectrum of a rotating magnetised star 

Armed with some idea of the mode spectrum of magnetised nonrotating stars, we next consider rotating magnetised config- 
urations. The earliest studies of magnetic oscillations suggested that the significance of the magnetic field on the oscillation 
spectrum should be linked to the ratio M/|W|; when additionally including rotational effects we would expect the relative 
significance of the two effects to be related to M/T. 

We first consider magnetic shifts in the /-mode frequency for rotating stars. Rotation splits the /-mode into co- and 
counter-rotating modes; we expect the frequencies of both branches of the mode to shift with the addition of magnetic fields. 
At low rotation, the magnetic shift for each piece of the /-mode is comparable with the shift in the nonrotating case, but at 
higher rotation rates the shift becomes less significant — see figure [5] This bears out our expectation that the magnetic shift 
should scale with M /T. 

We next turn to a-modes and r-modes of rotating magnetised stars. Based on our experience so far, we have expectations 
on how each mode should behave. We anticipate a rotational splitting of the a-modes into co- and counter-rotating pieces (as 
seen for the /-mode); in addition we expect to see some magnetic shift, scaling with M/T, in the r-mode. 

We begin by tracking the axial 2a-mode with increasing rotation, finding that as expected it undergoes rotational splitting 



Oscillations of rotating magnetised NSs with toroidal fields 11 




Figure 7. The m = 2,Zo = 3,4 hybrid magneto-inertial modes. Dashed lines represent the pure inertial (<B>= 0) modes, whilst solid 
lines show magneto-inertial modes, which reduce to pure Alfven modes in the Q — * limit. The left-hand plot shows the Iq = 4 (axial) 
hybrid modes, whilst the right-hand plot shows Zo = 3 (polar) modes. In each case the upper-frequency branch of a hybrid mode is seen 
to meet a corresponding i-mode as M/T — > 0. For the |ai mode, we were also able to track the lower-frequency branch, which appears 
to reduce to a zero-frequency mode in the M/T — > limit. 




Figure 8. In a slowly-rotating magnetised star, the r-mode is replaced by the axial Z = m a-mode. From the left-hand plot we see that 
this mode is not subject to the CFS instability if Q is sufficiently small, but at some higher rotational frequency /cfs ( a function of 
the field strength B) the a-mode crosses into the unstable regime. The right-hand plot shows the variation of /cfs with average field 
strength <B>. 



(figure[6]). The lower-frequency branch of this a-mode appears to tend to zero with increasing f2 (or equivalently, as M/T — > 0). 
The higher-frequency branch of the a-mode tends to the r-mode frequency as M/T — > 0. We confirm that the magnetic/inertial 
character of these hybrid modes depends on M/T by tracking the 2a-mode for three different field strengths, finding that 
when <B> is higher the hybrid-mode frequency approaches the r-mode frequency more slowly. The higher-frequency branch 
of the \a mode is counter-rotating and joins up with the r-mode, whilst the lower-frequency \a mode corotates with the star. 

Having established that the pure \a mode and the pure r-mode are replaced by a hybrid magneto-inertial mode when both 
magnetic and rotational restoring forces are present, one would expect to find similar hybrid modes corresponding to other 
Alfven/inertial modes; we confirm this expectation in figure [7] As before, rotation splits each a-mode into co- and counter- 
rotating branches. We are able to track the upper-frequency branches of both polar 2a-modes to their inertial counterparts, 
and all three |a-modes to known inertial modes in the M/T — > limit. In addition, we are able to track the lower-frequency 
branch of the |ai mode to high rotation rates; it appears to become a zero-frequency mode in the M/T — > limit, as for the 
lower 2 a- mode. 

A number of modes are subject to instabilities driven by gravitational radiation emission, but in general only become 
unstable for sufficiently rapid rotation. H owever, the r-modes are unstable even in very slowly rotating stars, in the absence 
of viscosity (jAndersson fc Kokkotasll200ll ). We have already seen that magnetic fields significantly alter the behaviour of the 
r-mode for slow rotation, so we now consider the effect this has on their stability. For a counter-rotating mode with frequency 
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Figure 9. A magnetic field changes the growth time of the r-mode instability. Here we plot an approximation of the ratio of magnetised 
{ t Gr)b to unmagnetised (tqr)o growth timescales, against dimcnsionlcss rotation rate. The dashed vertical line shows where the ratio 
asymptotes (i.e. when the magnetised mode becomes stable). We see that in all cases the instability growth is slower with magnetic 
effects, but the effect becomes insignificant for rapid rotation. The magnetic timescales used were for a star with a field of 2.87 X 10 16 G. 

a (positive by convention) in the rotating frame, the instability criterion is 

a(a - mft) < 0. (42) 

It follows immediately that radiative instabilities are suppressed when a > mfi. In the left-hand part of figure [8] we plot this 
threshold frequency, together with the nonmagnetic r-mode and the hybrid mode that replaces it in the magnetic case. It is 
clear that whilst the unmagnetised r-mode is always in the unstable regime, its magnetic equivalent (the hybrid of the r-mode 
and the axial I — m a-mode) is stable for sufficiently low rotation rates. The maximum rotational frequency a star can have 
before its \r mode goes unstable is plotted on the right-hand side, as a function of the stellar magnetic field strength. 

Even when magnetic fields are not strong enough to suppress the r-mode instability, they may slow down its growth. A 
full calculation of this effect is beyond the scope of this paper, but we may estimate it with some simplifying assumptions. 
The growth time tqr of the r-mode instability due to gravitational radiation is given by 

1 _}_dE 

tor 2E dt ( ' 

where E is the energy of the mode in the rotating frame. From this one can show that the growth time tqr scales with the 
rotating-frame mode frequency a in the following manner for an / = m r-mode: 



— a{a-m) 

tgr 



21 + 1 



(44) 



Andersson fc Kokkotas! (2001) for details. Note that for I = m = 2, the growth time scales with the sixth power of < 



We wish to estimate how the growth time for the ^r-mode changes when magnetic effects are included. Since tgw contains 
a factor of <r 6 , we will assume that this term has the most significant variation when a magnetic field is added. Other terms 
in the expression of E and its derivative will be approximated as constant. Using the indices and B to denote nonmagnetic 
and magnetic quantities (respectively), we then see that 

{tgr)b ^ oo(oo — 2Q)' 
{tgr)o ob(ob — 2Q) 5 

In figure [9] we plot this dimensionless quantity as a function of the rotation rate, finding that a toroidal magnetic field does 
indeed slow down the instability's growth. The importance of the effect depends on the rotation rate: at twice the threshold 
frequency for stability of the magnetised r-mode (i.e. when the mode is unstable), its growth time is still a factor of ~ 6 longer 
than in the nonmagnetic r-mode case; however, for very rapid rotation the difference in growth times is negligible. 



5 DISCUSSION AND CONCLUSIONS 

In this paper we have investigated oscillation modes of neutron stars with rotation and magnetic fields, specialising to the case 
of purely toroidal background fields; we intend to study purely poloidal and mixed-field geometries in future publications. Our 
numerical approach allows us to study oscillations of rapidly rotating and highly magnetised stars in a self-consistent manner. 
We first generate a stationary star in equilibrium to use as the background configuration; this star may have axisymmetric 
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distortions due to magnetic effects and rotation. We then time-evoive iinear perturbations on this background star in order 
to study its modes of oscillation. 

When a magnetic field is added to a star, the most obvious change to its oscillation spectrum is the presence of Alfven 
(a-) modes, a class of stellar oscillation restored by the Lorentz force. These modes are purely magnetic in nature only for a 
nonrotating background star. 

In a rotating magnetised star, we find that the pure a-modes of a nonrotating star (or equivalently, the purely inertial 
i-modes of an unmagnetised star) are replaced by hybrid magneto-inertial modes, w hose character is governed by the ratio 
of the magnetic M and kinetic T energies, as discussed by Morsink fc Rezania 1 20021 ). Tracking a star at fixed magnetic field 
from Q = through increasing rotation rate, we see a rotational splitting of the a-modes into co- and counter-rotating modes. 
The higher-frequency branches of these modes approach known i-mode frequencies. In general the lower- frequency branches 
are harder to track, owing to the dense nature of the oscillation spectrum, but when we are able to identify them we find that 
they appear to become zero- frequency modes in the M/T — > limit. 

The presence of these hybrid modes has parallels with other work. The evolutions of Passamonti et al. ( 2009T ) and 
Gaertig fc Kokkotasl (fcoosh found that when tracking p-modes (i.e. modes restored by composition gradients within the star) 



for increasingly rapid rotation, their frequencies approached known i-mode frequencies. One key difference between stratified 
and magnetised stars, however, is the behaviour of the r-mode in each case. Being purely axial in the slow-rotation limit, the 
r-modes are unaffected by composition gradients, whereas we have found that the presence of a magnetic field means that in 
the slow-rotation limit they become the axial / = m a-mode. 

Our work seems to be consistent with the analysis of iGlampedakis fc Anderssonl (|2007l ). who found that magnetic fields 
could act to suppress instabilities driven by gravitational radiation (the CFS instability) ; and in particular, that purely poloidal 
or purely toroidal fields should always play a stabilising role in this case. Using a to denote a mode frequency as measured in 
the rotating frame, it is known that modes satisfying the condition a [a — raQ.) < are susceptible to these radiation-driven 
instabilities; in particular, this includes the r-mode. In the presence of a magnetic field we find that the r-mode is replaced 
by the I = m axial a-mode; for sufficiently slow rotation we have a a > mfl and hence the mode is CFS-stable. In the regime 
where the star is unstable, we use a simple estimate to suggest that the instability's growth will be slower in the presence of 
a magnetic field. 

In addition to the hybrid magneto-inertial modes, there are also magnetic corrections to the /-mode frequency. These 
corrections are very modest (~1%) even up to field strengths of the order 10 17 gauss. In addition, as for the magneto-inertial 
modes, the magnetic correction becomes less significant still as M/T — > 0. However, we note that the magnetic correction 
seems to increase between m = 2 and m — 6, so although our approach limits us to low m one might expect more appreciable 
corrections to high-m /- and p-modes. 

Although it would be premature to make a quantitative comparison between our results and observed magnetar QPOs, 
we note that there are certain similarities in the oscillation spectra. The QPOs observed from SGR 1806-20 include 26 and 
30 Hz modes; these cannot be explained as overtones of crustal shear modes because the spacing is too small (they would 
need to be integer multiples for this). By contrast, it is easy to interpret these frequencies as global modes of a fluid star (i.e. 
the magnetar's interior), since we see modes at far smaller separation than integer multiples. For example, using our fitted 
relation (|4ip we see that the frequency ratio of the axial \a and polar \a\ modes is 0.030/0.033 ~ 0.91 — comparable with 
the observed ratio of 26/30 « 0.87. 

This pa per adds to the pic ture of magnetic stellar oscillat ions built up by a number of other recent numerical studies. 
The work of lsotani etafl (120081 ) and ICerda-Duran et al] (|2009l ) investigated axial magnetar oscillations, modelling the star's 
magnetic field as dipolar ( and hence pure l y polo idal). They found two localised families of QPOs, which they related to 



observed mag netar QPOs. IColaiuda etahl feoogh worked on a similar problem, but in the more general case of a mixed 



poloidal-toroidal background field. Their work complements o ther studies, but they are also able to identify a third family of 
QPOs in th eir model star. Finally, Sotani fc Kokkotasl ( 20091 ) find a set of polar oscillations of dipolar fields, agreeing with 
the work of iLed (|2007l ) that a magnetar should have both axial and polar oscillations. 

Many of these recent studies have analysed their results in the light of the suggestion that magneti c oscillations of a 
perfectl y-conducting sta r form a continuum, rather than di screte modes. This was proposed by Levinl d2007l). revisiting earlier 
work bv lGoossens et al.l (|1985l ). Various numerical studies (jSotani et alJl200Si : ICerda-Duran et a l. 2009; C olaiuda et al. 2009) 
have fo und results consistent with this proposal, in the case of axial oscillations of a dipole field. However, Sotani fc Kokkotasl 
( 20091 ) suggest that polar oscillations of a dipolar-field star are discrete. 

Since our background field is purely toroidal, we cannot make quantitative comparisons with work discu ssed in the last 



two p aragraphs, since those studies assumed dipolar fields (or mixed poloidal-toroidal fields in the case of IColaiuda et al 
(20091)). However, we do find broad similarities — in particular, our a-mode frequencies are of the order 100 Hz (for a field 
of ~ 10 16 gauss), as found from other magnetic evolutions. With a toroidal field there is no reason to expect a continuum 
of modes, since this feature has only be established for fields with a poloidal component. Indeed, all our results have shown 
discrete mode frequencies, with no dependence on position within the star, up to uncertainties due to resolution and the 
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finite duration of o ur simuiations (in practice, errors of ~ 
Sotani et al, (2008), but our axial a- modes are discrete too. 



1%). Our polar a- modes thus share this property with those of 



Purely toroidal fields and purely poloidal fields suffer from generic loc alised instab i lities, so in t he absence of damping 
mechanisms are not viable candidates for long-lived stellar mag netic fields <|Wrightlll97i lTavlerlll973h . Despite this, we have 
been able to perform stable evolutions of perturbations about a purely toroidal background for this work. There may be a 
number of reasons why these analytically-established instabilities have not affected our numerical work. Since we only consider 
first-order perturbations, higher-order effects are avoided; at the linear level, the greatest instabilities are those for m = 
and m = 1, whilst we have only considered m ^ 2. Finally, we have included artificial viscosity and resistivity to damp 
numerically-generated instabilities, and it is possible that these have prevented the growth of physical instabilities too . 

One way in which pure-poloid al/toroidal fields may be stabilised is through rotation ( Geppert fc Rheinhardt 20061 ; 
BraithwaitefeoOrj ; iKiuchi et al.ll2008 ). although this effect will be small in the case of the magnet ars, whose rotation al periods 



are very long. Relatively small poloidal components may stabilise dominantly toroidal fields ( Braithwaite! 200E), but it is 
difficult to draw general conclusions on the relative strengths of the two components, since other work ( Lander fc Jones 20091 ; 
Ciolfi et al. 2009T ) has found that apparently general constructions of magnetic stars in equilibrium (in both Newtonian and 



relativistic contexts) result in mixed fields which are dominantly poloidal. 

Given the many uncertainties regarding the nature of stellar magnetic fields, we believe that it is reasonable to study 
oscillations of purely toroidal fields, even though these may suffer certain instabilities, as we have discussed. Furthermore, 
a star whose field is dominantly toroidal could be expected to have an oscillation spectrum with similar features to those 
discussed in this paper. 

The observations of QPOs in the tails of giant flares from magnetars may be providing a rare probe of the interior fields 
of these stars. In addition, the longer-term prospect of studying gravitational- wave signals from neutron stars using detectors 
like Advanced LIGO and VIRGO should also help understand the physics of these stars. With a clearer picture of the nature 
of stellar oscillations in a strong magnetic field, we should be better equipped to interpret magnetar QPOs and GW signals 
of neutron stars. Our formalism allows us to build up a more sophisticated picture of the oscillation spectrum of magnetised 
stars by a gradual inclusion of extra effects. These could include other field geometries, stratified stars, a superfluid interior, 
and so on; we intend to explore these in future work. 
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